A numerical canonical transformation approach to quantum many body problems 
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We present a new approach for numerical solutions of ab initio quantum chemistry systems. The 
main idea of the approach, which we call canonical diagonalization, is to diagonalize directly the 
second quantized Hamiltonian by a sequence of numerical canonical transformations. 



INTRODUCTION 

The vast majority of current methods for numerically 
solving quantum many body systems fall into a few broad 
categories. One type of approach involves producing ap- 
proximate representations of wavefunctions. Examples 
of this approach include configuration interaction, cou- 
pled cluster methods [Q, and the density matrix renor- 
malization group (DMRG)|§, ||, |J. A second type of 
approach is based on perturbation theory. Some versions 
of perturbation theory are closely related to wavefunc- 
tion approaches, while others, for example utilizing finite 
temperature imaginary time Green's functions, are more 
closely related to path integrals. A third type of approach 
is quantum Monte Carlo, which also may involve repre- 
sentations of wavefunctions (e.g. Green's function Monte 
Carlo) or of path integrals (e.g. determinantal/auxiliary 
field methods). 

Much less explored than these approaches are meth- 
ods based on similarity or unitary transformations of the 
Hamiltonian H. In these approaches the primary focus 
is on H and transformed versions of it; wavefunctions 
play a much more minor role. In this paper we develop 
such an approach in the context of ab initio quantum 
chemistry calculations in a finite basis. This approach is 
based on unitary canonical transformations (CTs) of H 
written in second-quantized operator form. Such trans- 
formations have been used in analytical work for very 
long time|^]. A well known example in condensed mat- 
ter physics is the Schrieffer- Wolff transformation of the 
Anderson model into the Kondo model |6); another is the 
well known mapping of the Hubbard model in the large 
U/t limit into the t-J model. Often these transforma- 
tions are performed once, and relate one model system 
to another, simpler system, with fewer degrees of free- 
dom, in an approximate way. In some cases special CTs 
can produce exact solutions for certain model systems. 

Recently, substantial progress has been made by the 
development of continuous unitary transformations, in 
which a set of differential equations is solved to perform 
the CT. |?], |j| This method, which was developed indepen- 
dently by Wegner and by Glazek and Wilson, is known 
by the names "flow equation method" (?]] and "similarity 
renormalization" Q. A key advantage of this approach 
is that one does not need to know in advance the trans- 
formation operator to be used; it is determined implicitly 



by the solution of the differential equations. Another ad- 
vantage is that once the differential equations are set up, 
there is no operator algebra to be performed in the course 
of the numerical solution of the differential equations. 

This approach can be performed in a semianalytic con- 
text, where typically one is deriving one model from 
another. For example, an improved treatment of the 
Schrieffer- Wolff transformation of the Anderson model 
has been performed. |J It can also be used to obtain 
ground state energies and dispersion relations (i.e. exci- 
tation energies) for infinite lattice systems. 

Here we extend and develop the CT approach in a 
quantum chemical context. We first introduce a new 
version of the CT approach which is closely related to 
the Jacobi method of diagonalizing matrices. In this 
case, rather than solving a differential equation, one per- 
forms a sequence CTs, each involving the smallest pos- 
sible number of operators. This approach is designed to 
solve a finite quantum system, namely a molecule or clus- 
ter in a standard quantum chemical basis, as an alterna- 
tive to other standard approaches, such as configuration 
interaction or coupled cluster methods. We demonstrate 
that our method can determine ground states in ab initio 
chemical systems with excellent precision, utilizing tests 
on a water molecule for which exact results are avail- 
able. We expect similar performance for low lying excited 
states. We then present a version of the "flow equation" 
approach for use in the ab initio context. Utilizing ideas 
developed in our Jacobi method, we present a new version 
of the differential equations which make their numerical 
solution particularly efficient. We demonstrate that this 
method also works very well for the water molecule. 

Perhaps the most important aspect of our approaches 
is the ability to remove both low energy (i.e. core) and 
high energy virtual orbitals from the problem, leaving a 
system with a small number of "active" orbitals. This 
approach is especially useful for "strongly correlated" 
systems, i.e. those with open shells, breaking bonds, 
etc., where a single reference approach fails. In these 
systems, the strong correlation is generally confined to 
a relatively small number of orbitals. In wavefunction- 
based approaches, it can be awkward and expensive to 
deal with the strongly correlated part of the problem with 
a powerful method (e.g. full diagonalization of the ac- 
tive space) and the simple high-energy part of the prob- 
lem with another, simpler method. By working with the 
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Hamiltonian directly, one can separate the two parts of 
the problem simply and efficiently. First one transforms 
the whole Hamiltonian into a form in which the high and 
low energy orbitals arc either completely unoccupied or 
complete occupied (for core orbitals). Then these or- 
bitals can be thrown out of the problem completely, leav- 
ing a smaller system of partially occupied orbitals which 
can then be solved with any of a variety of nonpertur- 
bative techniques. Here we demonstrate this powerful 
hybrid approach on a stretched water molecule, solving 
the smaller systems with DMRG. 

Another advantage in working with the Hamiltonian 
directly is in obtaining excitations. In transforming the 
Hamiltonian, we make progress in solving for excited 
states at the same time as we obtain the ground state. Es- 
sentially the same approach is used to obtain both ground 
and excited states. Although we have not performed tests 
yet for excited states, CT approaches in other contexts 
have obtained excellent results. |I(J 

We will call the set of methods we present here canoni- 
cal diagonalization (CD). The term diagonalization indi- 
cates our intent to solve the system fully, rather than just 
transforming to a simpler model. The transformations in- 
volved are both unitary and canonical (which is defined 
in the next section); we choose only the term canonical to 
emphasize that the method works in the space of second 
quantized operators. We further distinguish the Jacobi 
CD method (JCD) and the flow equation CD method 
(FECD). All CD approaches share the feature that the 
object that one is manipulating is the second quantized 
Hamiltonian, as a collection of abstract operator terms 
with specific numerical coefficients. 

These approaches seem particularly suited to ab initio 
quantum chemistry, which are characterized by a very 
general quantum Hamiltonian, containing almost all pos- 
sible one and two-electron terms. The CTs generate ad- 
ditional terms involving one, two, three, and more parti- 
cles. Since general one and two particle terms are already 
present in the Hamiltonian, no extra inconvenience arises 
from these terms. Three and more particle terms are 
more inconvenient, but most such terms can be neglected, 
to an excellent approximation. CD is size-consistent, and 
many-particle terms which may be left out involve the 
simultaneous interaction of three or more (dressed) elec- 
trons, so that neglecting them is analogous to neglecting 
connected clusters involving triple and higher excitations 
in coupled cluster methods. Note that canonical trans- 
formations also appear in the theory of the coupled clus- 
ter method jyj, although the method remains largely a 
wavefunction approach. Note also that although one does 
not need to write any wavefunctions explicitly, CD in its 
simpler forms implicitly expresses the ground state us- 
ing the exponential of an operator acting on a reference 
state. Further links to coupled cluster methods are made 
in the Discussion Section. 

CD fits naturally into a renormalization group (RG) 



framework. First, one can remove ("integrate out") 
higher energy orbitals, one at a time if one wishes, leav- 
ing a system where the effects of the removed orbitals 
are incorporated into an effective Hamiltonian for the re- 
maining orbitals. Thus each step resembles a transforma- 
tion in a typical RG calculation in statistical mechanics, 
although unlike in that case one cannot continue indefi- 
nitely and there are no fixed points. Second, even if one 
is not integrating out orbitals, the transformation of the 
Hamiltonian, like in RG methods, occurs in a sequence of 
steps, with truncation of higher order terms occuring at 
each step. Third, the differential "flow" equation form of 
CD, in which a time- like variable controls the evolution of 
the Hamiltonian operator towards a more diagonal form, 
closely matches Wilson's original conception of the RG 
approach [^2). 

CD is a natural complement to DMRG, and this was 
a principle motivation in developing it. When applied 
to quantum chemistry problems, DMRG does very well 
in describing non-dynamical correlations (strong correla- 
tions associated with partially occupied orbitals), but it 
is inefficient in describing dynamical correlations, since it 
describes high energy virtual orbitals on the same foot- 
ing as partially occupied strongly interacting orbit alsfll3|. 
CD has complementary behavior. It can be used to re- 
move the nearly unoccupied orbitals, leaving a smaller 
Hamiltonian involving strongly interacting orbitals for 
DMRG to solve. 

JACOBI CD 

We begin with the Jacobi CD approach. In the ordi- 
nary Jacobi method for diagonalizing matrices, one ap- 
plies a large number of unitary transformations to a Her- 
mitian matrix to bring it into diagonal form|l4[. A uni- 
tary transformation gives a new matrix which has the 
same eigenvalues as the old. In the Jacobi method, each 
unitary transformations consists of a rotation of two rows 
and columns to zero out a single off-diagonal element Hij . 
The part of the unitary transformation matrix exp(^4) 
corresponding to rows and columns i, j 'is 

( e \ ( cose sme \ 

^ [-6 o) = {-sine cose) (1) 
The angle e which removes the term Hij is given by 

6=^taxr 1 [2H ij /{E i -E j )] (2) 

where Ei = Ha and the transformation is applied as 

exp(A)ff exp(-A). (3) 

Any unitary transformation can be written as an expo- 
nential of an antiHermitian matrix Ap5[|. For real sym- 
metric matrices and operators H and for what follows, 
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it suffices to use real antisymmetric matrices and opera- 
tors A. In the Jacobi method, one traverses the matrix 
repeatedly, rotating away off-diagonal elements, starting 
with the largest off-diagonal elements for efficiency. 

In Jacobi CD we construct unitary transformations to 
successively remove off-diagonal terms of a second quan- 
tized Hamiltonian. We consider a quantum chemical sys- 
tem in a Hartree Fock basis, with Hamiltonian 

H = 12 Ti i C L C 3<r + \ V ^klc\ a c\ a ,C ka ,C la . (4) 

ija ijklaa' 

Here T^- contains the electron kinetic energy and the 
Coulomb interaction between the electrons and the nu- 
clei, while Vijki describes the electron-electron Coulomb 
interaction. Indices such as i denote spatial Hartree- Fock 
orbitals, and a is a spin index. Later on we shall some- 
times use i to denote a spin-orbital, in which case the 
context should make the usage self-evident. We use a 
computer representation of H as a sum of abstract oper- 
ator terms, each with a coefficient and a product of c and 
c 1 operators involving various orbitals. In our program a 
c or operator is described by a single byte, which en- 
codes the orbital index involved, the spin, and whether 
it is c or c'. This implementation is thus limited to sys- 
tems with at most 64 orbitals. A complete operator term 
is stored as an array of such bytes plus a floating point 
coefficient. We developed a set of C++ routines to take 
products and commutators of such operators, putting the 
result in normal ordered form using the anticommutation 
relations. Having these operations be reasonably efficient 
is crucial to the method. This approach, using formal op- 
erator terms to describe H, rather than specific matrix 
expressions for terms of various orders, is both simple 
and general. In the future, in order to implement spe- 
cific approximations within CD, more efficient code could 
be produced by deriving and implementing the relevant 
matrix expressions. 

An off-diagonal term which can be rotated away is sim- 
ply any term which is distinct from its Hermitian conju- 
gate. Self-adjoint terms constitute the diagonal elements. 
They can always be written as products of density oper- 
ators ni a — c\ a Ci a . Consider, as a specific example, the 
term 

Vq = aV a = actj.ct^cfcjc;-]-, (5) 

where a is a numerical coefficient. Let 

A = 6{V a - V£) (6) 

Given the proper choice of 9, the transformation H — > 
exp(A)H exp(— A) will remove V a + from H, intro- 
ducing other terms instead. We will consider the choice 
of 9 momentarily. These additional terms in general have 
smaller coefficients than V a , making H more diagonal. If 
one continued the process indefinitely, one would have a 



"classical" Hamiltonian where every term was diagonal. 
In this case any Slater determinant ejet ... |0) is an eigen- 
state and all eigenvalues can be read off essentially by 
inspection. After the particle-hole transformation of the 
occupied orbitals described below, normally the ground 
state energy is simply the constant term in H . 

This unitary transformation is a canonical transforma- 
tion, in that when applied to the operators c and , the 
anticommutation relations are preserved, e.g. 

{e A Cl e- A , e A c]e- A } = e A {c u c]}e~ A = {c,, c}}. (7) 

One can view a CT as a complicated many-particle 
change of basis. In that sense CD is similar to the DMRG 
method, the biggest difference being that DMRG works 
in a wavefunction basis. In the Discussion section, we 
comment more on applying the CTs to the c operators. 

In order to carry out the transformations, it is conve- 
nient to use the well-known formula 

e A He- A = H+[A,H} + ±[A,[A,H]} 

+ [A, [A,H]]] + ... (8) 

A commutator between an TVi-particle term and an N%- 
particle term gives up to (N\ + N2 — I)-particle terms. 
Hence, if A is a one particle term, it does not generate any 
higher order terms. Indeed, rotations by one-particle A's 
correspond to single particle changes of basis. CD using 
one-particle A's can be used to perform a Hartee-Fock 
calculation, if the initial orbitals are not the HF orbitals. 
To evaluate Eq. (||), we generally treat each term in H 
separately. For most terms Vp the transformation has no 
effect, since [A, Vp] = 0. For relevant terms the expansion 
can be carried out order by order until all terms of a 
given order are neglible (i.e. their coefficients are very 
small). At each order the terms should be put in normal 
ordered form. The number of distinct terms generated 
from a single term Vp is a modest finite number, since a 
particular Ci or c\ can appear at most once in a term in 
its normal ordered form. 

For several reasons it is convenient to perform a par- 
ticle hole transformation on the occupied orbitals. Thus 
we define di„ — Ci a if orbital i is unoccupied in the HF 
state, and di a = c] CT if it is occupied. The d's have identi- 
cal anticommutation relations to the c's. After this trans- 
formation, anticommutation relations are used to put the 
terms in normal ordered form, putting d to the right of d' . 
The resulting anticommutators generate new lower order 
terms, after which the HF energy appears as a constant 
term in H. In terms of the d operators, the HF state is 
the vacuum |0). There are no off-diagonal single particle 
terms such as d\dj. There are terms which appear to vio- 
late particle conservation, such as d\d^d\d\, but which in 
fact do not in terms of real particles. As one performs the 
canonical transformations, the vacuum state approaches 
the exact ground state. 
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We now consider the choice of 9. An operator term V a 
connects an exponentially large number of pairs of states 
l,r together, (l\V a \r) ^ 0. However, one of these pairs 
of states can be considered the most important, namely, 
the one closest to the HF vacuum. This pair has the 
fewest number of rf^'s operating on |0) which generate 
states not destroyed by V a . As a specific example, let 
V a = 0.1 d\djdkd m . Then the most important pair of 
states is 

|r> = 44410) 

10 = 4|0>. (9) 

We will call these states simply the left and right states 
of V a . Other pairs of states, considered to be less im- 
portant, have additional dJ's, which do not appear in 
V a , applied to both |Z) and \r). For example, one could 
take the pair dJJZ) and d^r). Define for a term V a the 
left energy Ei = (l\H\l) and similarly the right energy 
E r = (r\H\r). We now choose to use the Jacobi formula 
for 9 to attempt to eliminate the off-diagonal term in 
the Hamiltonian connecting I and r. Thus, if a is the 
coefficient of V a , we choose 

e = ^tan- 1 [2a/(E l -E r )]. (10) 

This does not eliminate V a exactly, since the operator 
A connects many different states, not just I and r, so 
it is not a 2 x 2 transformation, as it would be for a 
matrix. Nevertheless, we find that this choice generally 
works well, typically reducing the size of the coefficient 
of V a by a few orders of magnitude. Note that the degen- 
erate case Ei — E r is nonsingular, generating an angle of 
±7r/4 (either angle can be chosen). Such a large trans- 
formation angle should be avoided if possible, however, 
since it generates high order terms in the transformation 
of H. 

A more common choice in analytic work using CTs is 
to choose to eliminate V a to first order in the expan- 
sion Eq. (U), namely choosing 9 to set the coefficient 
of V a in [A, H D ] to —a, where H D is the diagonal part 
of the Hamiltonian. This is closely related to our ap- 
proach: note that (l\ [V a , H D ] \r) is the coefficient of V a 
in [V a ,H D ]. However, 

(l\[V ai H D ]\r) = (E r - Ei)(l\V a \r) = E r - E l (11) 

so that this choice gives 9(E r — Ei) = —a. This agrees 
with our choice to lowest order, but it is not well behaved 
if E r w E t . 

One need not eliminate all off-diagonal terms in H. 
If one is interested in only the ground state, then one 
needs to eliminate all terms connecting that state to other 
states. More specifically, suppose the initial HF state 
|0) has substantial overlap with the ground state. Only 
states which produce a nonzero result when acting on 



the vacuum need be removed, namely, only terms such 
as djd k dl n dl l or djd k , and their Hermitian conjugates (as 
well as similar multiparticle terms). Since all terms still 
satisfy particle conservation, in each of the two-particle 
terms two of the orbital indices jkmn must correspond to 
occupied states, and two to unoccupied orbitals. Once all 
such terms are eliminated, then the state |0) is the ground 
state, and the ground state energy is the constant term 
in the Hamiltonian. 

In all but the smallest systems, some of the terms 
formed from the CTs must be discarded according to 
some criterion. The simplest criterion is to neglect all 
terms involving three or more particles, i.e. six or more 
d operators. Our test calculations on the water molecule 
suggest that this is a very accurate approach for systems 
well described by a single reference state. Other possible 
criteria include keeping all terms whose coefficients are 
larger than some cutoff; keeping all one and two particle 
terms and all three particle terms larger than a cutoff, 
etc. More sophisticated criteria are possible also, such 
as trying to estimate the contribution of each term using 
pertubation theory, and discarding terms whose contri- 
bution is below a cutoff. Here, we perform some test 
calculations according to simple cutoff criteria. In the 
future, we hope the criteria can be optimized. 

In order to preserve symmetries, such as spin symme- 
try, one can rotate sets of terms which are related by a 
symmetry transformation in one step. For example, in 
what follows, for each term, we check to see if it is dis- 
tinct from the term coming from flipping all of its spin 
indices. If it is distinct, both are rotated together with 
the same rotation angle. The rotation angle for both is 
chosen as the angle to rotate one of the terms separately. 
This procedure preserves spin symmetry exactly. 

One must decide in which order to go through the 
terms in performing the CTs. Since each CT alters the 
coefficients of many other other terms in H, it makes 
sense to start with the largest first. One approach would 
be to find the term with the largest magnitude coefficient 
at each step. Another would be to choose the largest ro- 
tation angle. However, searching for the largest term 
at each step would be inefficient. Therefore, we have 
chosen the following method: a cutoff angle is chosen, 
and all terms with angles greater in magnitude than this 
cutoff are treated in a sweep through the terms, in a 
predetermined but arbitrary order. Then, the cutoff an- 
gle is reduced by a constant factor, and the procedure 
is repeated. Here, we started with a cutoff of 0.15 and 
reduced it using a factor of 0.6. (In some passes, partic- 
ularly the initial one, there may be no CTs performed.) 
In Table I we show results for a 25 orbital DZP basis wa- 
ter molecule, for which full CI results are available Jl6[. 
Because of some arbitrary choices in the ordering of the 
CTs, which unfortunately can affect the results slightly, 
the results here would be difficult to reproduce precisely 
by an independently written program. (The differen- 
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TABLE I: Results from the Jacobi CD method applied to 
the water molecule in a 25 orbital basis. The Is O core HF 
orbital has been frozen. The exact energy of the system in 
this basis is -76.256624 Hartrees. In all cases, all two particle 
terms have been retained. £3 and £4 are the cutoffs for re- 
taining three and four particle terms, and -/V3 and Na are the 
corresponding maximum number of such terms in H during 
the diagonalization. AE is the error in the energy, E — i? e xact- 
AE N 3 N 4 



00 00 

0.01 00 

0.001 00 

0.0005 00 

0.0001 00 

0.0001 0.0001 



0.0041 

0.0041 

0.0019 

0.0011 

-0.0001 

-0.0001 





464 

1.1 x 10 5 
3.6 x 10 5 
2.9 x 10 6 



2.9 x 10 6 2.0 x 10 s 



-76.0 



-76.1 



-76.2 




-76.3 



• • Constant Term 

* x Pert. Thy. 

Exact 



»««»»» • • _ 



500 1000 1500 

n 



tial equation method discussed below does not suffer this 
problem.) Despite this, one can easily evaluate the po- 
tential of the method from our results. One can see that 
we obtain an accuracy of several millihartrees even when 
truncating all three or more particle terms. If one keeps 
also some three particle terms, one can obtain accuracy 
to fractions of a millihartree. This sort of accuracy is 
comparable to coupled cluster methods. 

It is not necessary to perform every CT in this proce- 
dure. At the end of every sweep, we perform a calcula- 
tion of the energy using second order perturbation the- 
ory for the current Hamiltonian. The calculation time for 
this procedure scales only as the number of terms in H, 
so there is a neglible impact on the overall computation 
time. As the largest angle terms are eliminated, the per- 
turbation theory result becomes more and more accurate. 
Even just a few rotations of the largest terms can make 
perturbation theory much more well-behaved. One can 
stop the procedure when the perturbation result is well 
converged, which typically happens long before all the 
chosen terms are removed. In Fig. 1 we show the results 
for this procedure. One can see that the perturbation 
result converges much more quickly than the constant 
term in the energy. One might well stop after about 300 
Jacobi steps; in this particular example this number is of 
order N 2 , where N is the number of orbitals. 

In Table II we show similar results for a water molecule 
whose bonds have been stretched by a factor of two. This 
system is not well described by a single reference state: 
in the full CI calculations of Olsen, et. al.|l7j|, in a differ- 
ent but similar basis, the weight of the HF determinant 
in the full CI wavefunction was 0.589, versus 0.941 for 
the unstretched molecule. Here, we find that CD is un- 
stable if only two-particle terms are kept. One finds that 
repeated Jacobi diagonalization steps reduce the energy 
without bound. CD is exact if no truncations are made, 
so this is an artifact of the truncation of three and more 
particle terms. However, keeping even a large number 
of three particle terms does not result in a particularly 



FIG. 1: Energy for the water molecule of Table 1 as a func- 
tion of the number of Jacobi rotations performed n. At each 
sweep the constant term of H is shown, as well as the cur- 
rent result from second order perturbation theory. The initial 
value of the constant term is the HF energy; the initial value 
of the perturbation theory is what one would get from it with- 
out doing CD. 



TABLE II: Same as for Table I, but for the water molecule 
with bond stretched by a factor of two. The exact result is 
-75.95227. 

£3 £4 AE A 3 A 4 
00 00 —00 
0.01 00 0.007 1.6 x 10 4 
0.001 00 0.015 1.9 x 10 5 
0.0005 00 0.015 3.8 x 10 5 
0.01 0.01 0.032 1.3 x 10 4 3.9 x 10 3 



accurate calculation. 

For this system, examination of the occupancies of the 
HF orbitals in the exact ground state (which we have 
computed with high accuracy with DMRG) reveals that 
there are four spatial orbitals with occupancies far from 
or 2; specifically, they have occupancies of 1.58, 1.52, 
0.46, and 0.4. The rest have occupancies less than 0.03 
or more than 1.97. In the case of the unstretched wa- 
ter molecule, occupancies are all within 0.05 of or 
2. The results for occupancies of natural orbitals are 
very similar |L7j. The contribution to the energy of a 
Hamiltonian term A, (ip\A\ip), can be expressed as a 
Green's function or density matrix element. In the case 
of Hamiltonian terms made out of operators involving 
only nearly filled or unfilled orbitals, the behavior of the 
Green's function is well understood, and the magnitude 
falls rapidly as one considers terms involving more parti- 
cles. For partially occupied orbitals there is no reason to 
believe that three or more particle Green's function ele- 
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TABLE III: Same as for Table II, but for a truncation crite- 
rion with no limit on the number of partially occupied terms. 
Here A3-1- is the maximum total number of terms involving 
three or more particles. Here the Is orbital has not been 
frozen; the "exact" value is taken from a DMRG calculation 
keeping 750 states: -75.9661. 



AE N: 



3+ 



0.001 0.0097 5.3 x 10 5 



0.0005 0.0025 1.4 x 10" 
0.0002 0.0003 3.9 x 10' 



,6 



merits are small. Consequently, one should only truncate 
such a term if its coefficient is small. For this reason, 
we have performed test calculations with the following 
truncation criterion: all terms with more than four d and 
d' operators corresponding to non-partially-filled orbitals 
are truncated. In addition, all terms whose coefficient's 
magnitude is below a cutoff e are eliminated. If there are 
N p partially filled orbitals, then this rule allows terms 
with up to AN p + 4 d's to appear. In this case we have 
up to 20 d's, i.e. a 10-particle term. 

As shown in Table II, with this criterion we see sub- 
stantially better results: we find that in this non-single 
reference system, accuracy to fractions of a millihartree 
is possible. 



FLOW EQUATION CD 

It is also possible to formulate CD in terms of a dif- 
ferential equation. This approach was originally devel- 
oped independently by WegnerQ and by Glazek and 
Wilson, Q in rather different contexts than we present 
it here. We will derive it here as a natural variation of 
the Jacobi CD method. In this approach we introduce 
a time-like variable i, and the Hamiltonian evolves as t 
increases. First consider a time-dependent Hamiltonian 
for some fixed antihcrmitian operator A: 



Hit) 



v H(0)e 



Here H(0) is the initial HF Hamiltonian. We have 
dH(t) 



dt 



= [A,H(t)]. 



(12) 



(13) 



This differential equation form of a CT has long been 
used in analytical work, where one integrates t from 
to some fixed rotation angle. Here we modify this by 
making A depend on H. First expand H as follows 



H(t) = ^a a (t)h a . 



(14) 



Each h a is a product of d and d' operators, and a a (t) is 
the corresponding coefficient. Let 



A(t) = 22 s a a a (t)h a . 



(15) 



The s a are fixed parameters, which we initially consider 
to have only three possible values: ±1, and 0. We set s a 
to if we are not interested in rotating the coefficient of 
h a to 0, because, for example, h a does not act directly 
on the HF state |0). For terms we wish to rotate to zero, 
we choose the sign of s a so that (1) A(t) is antihermitian, 
and (2) increasing t rotates in the direction to diminish 
a a (t) . These conditions are satisfied if s a is chosen as the 
sign of Ei—E r , where I and r are the left and right states 
of h a . We evolve H (t) as a sequence of infinitesimal CTs, 
as follows 



H{t + 5t) = e StA ^H{t)e- StA ^ 

= H(t)+St[A{t),H{t)}+0{St) 2 . 



(16) 



In the limit that 5t — > 0, this is equivalent to solving the 
nonlinear differential equation 



dH{t) 
dt 



[A(t),H(t)}. 



(17) 



Each infinitesimal rotation acts to diminish each a a (t) 
with nonzero s a . Since Ait) depends linearly on the 
a a (t), the rotations become smaller as the a a {t) decrease. 
Thus, we expect the solution of this equation for t — > oo 
to have a a (t) = if s a is nonzero. We also expect these 
a a {t) to diminish exponentially with t. 

If no truncations are made, the solution to this dif- 
ferential equation for any time t gives an H(t) which is 
related to H by an exact CT. This is true also for any 
choice of the s a as long as they satisfy the requirement 
that if hp — h' a , then sp — —s a , ensuring that A(t) is 
antihermitian. For numerical efficiency, it is useful to 
modify the choice of s a . This is because different terms 
h a require different rotation angles. One would like to 
make the exponential decay to zero of each a a (t) have ap- 
proximately the same time constant. If they have widely 
varying time constants, the number of steps in integrat- 
ing the differential equation will be very large. We can 
achieve this by choosing, for nonzero s a , 



Sq — (Ei — E r 



(18) 



Provided a a (t) <C Ei — E r , this choice makes the coef- 
ficient of h a in A(t) the angle 9 required to rotate the 
term to zero. This makes the natural time scale for each 
term equal to unity. We choose the s a at the beginning, 
using the untransformcd HF energies, and never change 
them; however, one could also make the s a depend on t. 

In Wegner's original flow equation method, rather than 
the above forms of A defined in terms of s Q , one took 
A = [H D , H], where H D is the diagonal part of H. This 
is very similar to the choice s a = Ei—E ri assuming all off- 
diagonal terms are being removed. However, this choice 
gives very widely varying time scales, driving terms with 
large Ei — E r to zero much more quickly. In the sense 
that the large energy difference terms are removed first, 
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Wegner's method can be considered a renormalization 
group method in itself, and one might stop at some finite 
time and study the partially transformed Hamiltonian. 
Our choice is much more efficient numerically, assuming 
one only wants the t — > oo limit. 

We decribe all the commutator relations in terms of a 
"matrix" B 



(19) 



If a commutator gives a term which is not in the set 
of Hamiltonian terms we are keeping, then that term is 
ignored. Then the final form for the flow equation CD 
method is a set of differential equations 



da-y(t) 
dt 



(20) 



a/3 



which are to be solved numerically. The B matrix was 
computed initially and stored in our program. Because 
of some regularities in the pattern of nonzero elements 
of B, the storage could be reduce by a factor of about 
N, the number of orbitals, from a naive estimate. How- 
ever, they could also be recomputed at each step to save 
storage, at the expense of computer time. Another ap- 
proach to save storage would be to remove a few orbitals 
at a time. One could even remove one term at a time by 
making only one s a nonzero, in which case the flow equa- 
tion method becomes very similar to the Jacobi method. 
To integrate the coupled differential equations, we use 
a simple fourth order Runge Kutta method with auto- 
matic step size adjustment. This routine attempted to 
integrate the differential equations with an absolute error 
tolerance of 10 -8 , and we integrated the equations from 
t = to t = 20. 

In Figure ||, we show the evolution of the constant 
term in H as a function of t for the unstretched water 
molecule. Only one and two particle terms were retained. 
The step sizes used were rather large, and they steadily 
increased. They are visible via the circles in the curve. 
Only twelve steps were taken, although each RG step 
in our very crude integrator required twelve derivative 
evaluations, Eq. (^p|). The result for the energy was in 
error only by about a milli-hartree. 

In Figure ||, we show similar results for the stretched 
water molecule. As in the Jacobi method, CD keeping 
only two particle terms is unstable, with the energy tend- 
ing to — oo. We believe that by keeping multiparticle 
terms one could make this method perform very well on 
the stretched water molecule, just as we found for the 
Jacobi method. 



REMOVING SETS OF ORBITALS 

Another approach for systems such as the stretched 
water molecule, which have some strongly correlated or- 



-76.0 



-76.1 



ha 




-76.2 



-76.3 



FIG. 2: Evolution of the constant term in the renormalized 
Hamiltonian as a function of time, for the flow equation CD 
method. The system is the same as in Table I. All two particle 
terms were retained in H. The final energy is -76.25795, versus 
the exact full CI value of -76.25662, shown by the dashed line, 
for an error of 1.3 milli-hartree. 




FIG. 3: Same as for Fig. g, but for the stretched water 
molecule. Here, the flow CD method retaining only one and 
two particle terms is unstable. 



bitals, is to integrate out many of the non-strongly- 
correlated-orbitals, leaving a small but strongly corre- 
lated system to solve with CD retaining many-particle 
terms, with DMRG, or with another method. We first di- 
vide the orbitals into two sets, those to be kept and those 
to be removed. Some of the orbitals to removed will have 
occupancies near 0, and some may be core orbitals with 
occupancies near 2. Due to the particle- hole transforma- 
tion, we need make no distinction between these cases. 
Consider the many-particle basis states \s) = d\d\ . . . |0). 
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Let \s) denote all states in which no orbital to be removed 
is occupied. Conversely, let \s') denote the rest, in which 
at least one orbital to be removed is occupied. We wish to 
rotate away all Hamiltonian terms which connect states 
\s) to \s'). Let r represent orbitals to be removed. Then 
the terms to be removed are described by the following 
rule: the terms have one or more dps, or one or more 
rf r 's, but do not have both dps and rf r 's. 

One finds that typically a few of these terms to be re- 
moved, largely by accident, have E r and Ei nearly iden- 
tical, although neither is close to zero because they in- 
clude operators adding or removing high-energy orbitals. 
To remove these problem terms requires a large angle of 
rotation. This can be disastrous for either the Jacobi 
or flow equation method unless many-particle terms are 
kept. However, some reflection indicates that these prob- 
lem terms are likely to be quite unimportant in terms of 
their true contribution to the ground state. Consider an 
nth order perturbation theory contribution to the ground 
state energy. Ignoring the energy denominators, such a 
term is proportional to 

(0\h 1 h 2 ...h n \0). (21) 

Clearly, if this contribution is nonzero, term h n must have 
E r = 0, and term h\ must have E\ = 0. Thus our prob- 
lem term with nearly degenerate nonzero energies cannot 
contribute in second order. For third order, one might 
consider hi to be the problem term. However, for this 
term either |r) or |Z) must belong to the set \s'). Remov- 
ing all the other terms connecting \s) to \s') means that 
either hi or /13 is to be removed, since |0) belongs to |s). 
Thus there is no third order contribution. The lowest 
order contribution for such a term is in fourth order, in- 
volving both this term and its Hermitian conjugate as hi 
and or involving two such terms as hi and ft.3. Here 
hi takes one from |0) into a higher energy state |s), /13 
takes one from \s) to |s'), hi then takes one back to |s), 
and hi takes one back to |0). The energy denominators 
are well-behaved, since E r and Ei of the problem term 
are not close to zero. 

Thus a quite reasonable approach is to rotate away all 
terms connecting \s) to |s') except those whose energy 
difference \E r — Ei\ is below some cutoff d. The terms 
below the cutoff are retained during the CD process, dur- 
ing which time they may change due to other terms being 
rotated away. After the CD is complete, one then can dis- 
card all terms having any d r or d\ operators, which will 
include these problem terms. One can also solve the CD- 
transformed H before truncation, using another method, 
and check that the occupancies of the removed orbitals 
are very close to zero. 

In Table IV we show the results of such calculations for 
the stretched water molecule. There are three sources of 
error in these calculations. First is the DMRG error, typ- 
ically near 0.0002 mH, keeping 400-600 states, which is 
small enough to show the other sources of error. Second 



TABLE IV: Results for the flow equation method applied to 
integrate out a set of orbitals, coupled with DMRG to solve 
the resulting Hamiltonian. The system is the stretched water 
molecule of Table II, with 25 orbitals. The first column tells 
how many orbitals were removed, specified as having the high- 
est single particle energies in the particle-hole-transformed 
Hamiltonian. The parameter d is the lower limit on the en- 
ergy difference of an operator for it to be removed. A_Eqd is 
the error in the energy, as computed by DMRG, relative to 
the full CI energy (-75.95227) after CD has been performed 
to eliminate orbitals, but with all 25 orbitals still present. 
n mal is the largest occupancy of any of the orbitals which 
have been "rotated away". A_Bcdt is the error in the en- 
ergy, computed by DMRG, after CD and after truncation of 
the rotated orbitals. The * indicates that the ground state in 
this diagonalization has a clearly erroneous orbital occupancy 
pattern, indicating that it is a low lying excited state which 
has dropped below the true ground state. The true ground 
state occupancy pattern reappeared upon truncation of the 
rotated orbitals. 
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0.011 
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0.012 



is the error from performing CD keeping only one and two 
particle terms. This is given by AEqd. Increasing d, or 
removing fewer orbitals improves AEqd. Third, there is 
the energy from throwing away the removed orbitals after 
CD. This is measured by the difference between AEcn 
and A£?cdt, and also by the maximum occupancy of 
the removed orbitals n max . We find that d can be made 
quite large: 0.5 is always fine, whereas 1.0 can be too 
large. We also find that we can remove up to about one 
half of the orbitals and incur only a very small error, 
even only keeping one and two particle terms. For the 
resulting small system even full CI would be a very easy 
calculation. Even removing all but four of the orbitals 
we get a reasonable result. We have not carried out any 
similar calculations keeping many-particle terms, but we 
can deduce the probable outcome. Since all the rotation 
angles 9 are rather small in this procedure, four particle 
terms, which can come in only as 9 2 , would be negliblc. 
Three particle terms come in as 9, and if such a term 
only involved the retained orbitals it presumably would 
have both Ei and E r small and it could give a substan- 
tial contribution to the energy of order 9. Three particle 
terms involving removed orbitals would have Ei or E r 
reasonably large, and would only contribute to the en- 
ergy via second order perturbation terms, thus coming 
in as 9 2 , which could be neglected. In short, we expect 
that keeping three particle terms involving the retained 
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orbitals only would be a very accurate approach for re- 
moving more than half of the orbitals. 

We would like to conclude this section with an argu- 
ment that the proper way to separate the treatment of 
high-energy from low-energy orbitals is by using an effec- 
tive Hamiltonian to remove the high energy orbitals, as 
we have done, rather than any wavefunction based ap- 
proach. We will make this argument via a trivial 3x3 
matrix, designed to have some of the crucial features of 
a strongly correlated/multireference system. Define the 
matrix 



H(e,S) 




(22) 



The third row and column represent a high energy or- 
bital, which we would like to treat separately from the 
first two nearly degenerate rows and columns. We will 
consider the parameter values (e = 0.1, S = 1), (e = 0.1, 
5 = 0.5), and (e = 0, S = 1). For these three param- 
eters we find the following ground state energies and 
eigenvectors (respectively): -0.099, and (0.995,-0.0098,- 
0.098); -0.064, and (0.789,-0.614,-0.022); and -0.196, and 
(0.700,0.700,-0.137). Now suppose we wanted to solve 
this system in two steps, first treating the third "orbital" , 
then next the other two, using a wavefunction approach. 
In treating the third orbital we insist that we ignore the 
small parameter e; otherwise we are treating the whole 
matrix together. We imagine that we have some pertur- 
bative method for obtaining the third component of the 
wavefunction, ignoring e; with this fixed, then we obtain 
the first two components, taking e into account. How- 
ever, comparing the first and third sets of parameters, 
we see that the third component ip 3 depends strongly on 
e, so this method must fail. 

Alternatively, we might imagine first treating the first 
two rows and columns separately, ignoring S and finding 
the ratio of components ipx/ip2> and then subsequently 
using 6 to fix ip 3 . In this case, comparing the first 
and second parameter sets, we see that ipi/^ depends 
strongly on 5, so that this method fails. In short, to treat 
this problem successfully, wavefunction based approaches 
must treat both S and e simultaneously. 

Now consider a simple CT approach. Rather than us- 
ing the Jacobi or flow equation method, we use a less 
sophisticated, but well-known perturbative CT method 
for removing the third row and column. Jl^| In this case, 
we find that the second-order change in the upper left 
2x2 portion of the matrix, due to the third row and 
column, is 



AHtj - H l3 H 3j -(- 



1 



2 Ei — E 3 



Ej - E 3 



) 



(23) 



removed k.) e appears only in the energy denominators, 
as a small correction; we ignore it by setting it to zero 
there. We obtain H cS = H + AH as 



H cS {e,5) 




(24) 



where Ei = Ha. (The general formula is obtained by 
replacing 3 by k and summing over all orbitals to be 



The ground state energies and eigenvalues for H eS for 
the three cases are -0.1, and (1.0,0.0); -0.064, and (0.788,- 
0.615); and -0.2, and (0.707,0.707). These results com- 
pare very nicely to the exact results for the full matrix. 
Indeed, they must; the procedure is well controlled, with 
large energy denominators. 

In order to properly separate the two parts of the prob- 
lem in a wavefunction-based approach, one needs to allow 
a set of possible wavefunctions to represent the high en- 
ergy states, rather than a single part of a wavefunction. 
Such an approach is embodied in the DMRG method, 
which chooses the optimal set of states to represent each 
part of the system. 



DISCUSSION 

CD is size-consistent: if one duplicated the Hamilto- 
nian for a system, corresponding to having two molecules 
separated by a large distance, and put in no interaction 
terms between the two systems, then no interaction terms 
would ever be generated and each system would behave 
identically under the CTs. The energy would be double 
the energy for one system. 

The calculation time for CD generally scales identi- 
cally with the number of orbitals N for the Jacobi and 
flow equation methods. Consider first the methods which 
directly determine the ground state, rather than remov- 
ing orbitals first. There are of order N^ cc N^ nocc terms 
r which connect directly to |0), where N occ (N^ nocc ) are 
the number of occupied (unoccupied) orbitals, which one 
needs to remove. Not all other terms s connect to any 
term r; if one is discarding all three particle terms, then 
there must be two orbital indices matching in r and s to 
get a contribution. Thus each term r connects to of order 
./V 2 terms s. Hence the total calculation time scales as 
N^ cc N^ nocc N 2 , or roughly N^ CC N 4 or more roughly N 6 . 
This is comparable to a singles and doubles CI or coupled 
cluster. If one treats only the terms r with large angles, 
using second order perturbation theory for the rest, the 
calculation time would be reduced but the scaling is more 
difficult to analyze. However, from the results of Fig. 1 
it is tempting to estimate the number of terms needed 
to be rotated as about iV 2 , leading to an overall scaling 
of TV 4 (plus a time of order N 5 for the initial HF change 
of basis). Of course, studies of systems of various sizes 
are necessary to determine the true dependence on N. 
(It is also challenging to write efficient programs for CD 
which exploit the potentially favorable scaling: if one is 
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not careful, one may find one's program spending most 
if its time performing commutators very slowly for terms 
with small coefficients which are later discarded.) One 
could also only rotate the largest N 2 terms using the 
flow equation method, and then use perturbation theory 
for the rest of the terms, leading to similar scaling with 
system size. There are also other variations of CD with 
good scaling. Note that if one does CD but restricts the 
terms s to be either r ' , or a diagonal term whose indices 
all match those in r, then one obtains an o(N 4 ) method 
closely related to second order perturbation theory. A 
presumably more accurate o(N 5 ) method is obtained if 
one restricts s so that three indices must match those in 
r, rather than two. For CD where one removes sets of or- 
bitals, keeping one and two particle terms, the scaling to 
remove each orbital is o(N 2 cc N 3 ), for a total of N 2 CC N A 
to remove a finite fraction of the orbitals. 

Let us discuss in more detail how to think about the 
canonical transformations jH^. Thus far, we have taken 
the view that we apply a CT to get a new Hamiltonian 



If we define new operators d using the inverse CT, 



H = e A He 



(25) 



which has different coefficients from H, but is written in 
terms of the same operators 



H 



(26) 



The new Hamiltonian has the same eigenvalues as the 
old, and one can reconstruct the eigenvectors: if 



then 



H\1>) = E\1>), 
He- A \ij))=Ee- A \il;), 



(27) 



(28) 



so that e A \ip) is the corresponding eigenvector of H. 
One could also define new operators di and d\ as 

di = e A d t e- A , (29) 

where the same expression applies for d\. Since 

e A didje~ A = e A d i e~ A e A d :j e~ A = didj, (30) 



one could equally well write H as 
H = 2^ a aha- 



(31) 



Here h a is a product of di operators with the same orbital 
indices and order as h a . 

This form, Eq. (31), is not especially useful, since the 



coefficients of the Hamiltonian are not any more diagonal 
than in H. A more useful expression comes from writing 



d t = e- A d t e A , (33) 



then 



(34) 



where h a is defined analogously to h a . We see that in 
terms of the d operators, the original Hamiltonian has the 
more diagonal form for the coefficients of H. This means 
that one should think of d\ as the operator which creates 
a quasiparticle, not d\. In particular, suppose A fully 
diagonalizes H, in which case any Slater determinant is 
an eigenstate of H. For any orbital i 



Hd\\Q)=e i d\\V), 



from which we obtain 



Hdle~ A \0) =E t dle~ A \0). 



(35) 



(36) 



H 



- A e A He~ A e A 



'He' 



(32) 



We see that d\ creates a new exact eigenstate from the 
ground state e~ A \0), containing an extra particle associ- 
ated with orbital i. This defines df to be a quasiparticle 
creation operator. It creates a "dressed" electron, with 
correlations built in. Because of the correlations built in, 
three and more particle terms can appear in H. Note 
that if one has exactly diagonalized H with A, then one 
can create all of the excited states by successively apply- 
ing d\'s to e~ A \0). 

The formulation of CD in terms of exponentials of op- 
erators has much in common with coupled cluster meth- 
ods (CC). In coupled cluster methods, the ground state 
wavefunction is written as e T |0). Usually T is not anti- 
hermitian, but in some less common versions of CC, it 
is, and usually the CC equations are derived using (for- 
mally) a similarity transformation of H.jlJ One difference 
between the two is that in CD we never explicitly write 
down A] rather, we perform a sequence of transforma- 
tions A%, A2, ■ ■ ■ A n , which implicitly define the complete 
transformation e A = e An . . . e Al . (In the flow equation 
method this sequence is continuous.) Based on the simi- 
lar expressions for the ground state, one might expect CD 
and CC to have similar errors, and our results are gen- 
erally consistent with this. However, the overall point of 
view between CD and CC is fundamentally different: CC 
is approximating the ground state, whereas CD is pro- 
gressively transforming the Hamiltonian into a diagonal 
form. The point of view of CD makes certain approaches 
natural and manageable, including removing sets of or- 
bitals, extracting excited states, and utilizing renormal- 
ization group ideas. 

Furthermore, CD, in its various approximate forms, 
makes its truncations of H at each transformation. These 
intermediate truncations make tractable the use of uni- 
tary transformations, rather than non-unitary similarity 



11 



transformations. Such continuous truncations are famil- 
iar from RG methods in statistical mechanics. One way 
of understanding their usefulness is to consider diagonal- 
izing a matrix with an approximate second order uni- 
tary transformation, as in the previous section. Here, 
however, we consider transforming the whole matrix this 
way. Except for matrices which are nearly diagonal to 
start with, this second order approach would work very 
poorly. However, if one makes a sequence of second-order 
unitary transformations, each having very small rotation 
angles, the method becomes accurate; in fact, it is exact 
in the continuous limit. This is analogous to integrat- 
ing an ordinary differential equation very precisely with 
a sequence of very small time steps, using a low order 
integration method. This is also why the flow equation 
CD method, without truncation, is exact even though 
only a first order commutator appears in the equation. 
The truncation of many particle terms is not really anal- 
ogous to throwing away higher order commutators, and 
so CD with truncation is not exact. However, there is 
no reason a priori to expect that CD, with its continous 
truncations, should be worse than CC. 

Now let us briefly mention how to obtain excited states. 
Suppose one wants to know the energy of an excited state 
which has a large overlap with the state d\\0). One needs 
to remove all off-diagonal terms which do not destroy 
this state, such as d^d^djdi, plus their Hermitian conju- 
gates. This includes terms such as d^d^jd^, which one 
would already remove to get the ground state. It may 
happen that some of these new terms to remove would 
require large rotation angles, in which case one might 
want to remove most of the orbitals first. Note that if 
one removes a large number of orbitals, a full diagonal- 
ization obtaining all excited states of the remaining or- 
bitals may be quite manageable. One might also try to 
remove all off-diagonal terms in H , in which case all the 
excited state energies could be obtained by inspection! 
Note that the work for removing all off-diagonal terms in 
H would still scale as N 6 . However, in this case, there 
would be many terms with nearly degenerate E r and E[ 
which would cause problems. We leave exploration of 
these approaches for future work. 

Let us also briefly mention calculation of expectation 
values of operators in the ground state, (A). One ap- 
proach is simply to apply the same CTs to A as one has 
applied to iJ, truncating many-particle terms in a simi- 
lar fashion, to get A, and then evaluate (0|A|0). Another 
approach would be to obtain an approximate expression 
for the ground state in the original HF basis, by ap- 
plying cxp(A) successively to |0) for each CT in reverse 
order, again with some truncation rules. Again, we leave 
exploration of these approaches for future work. 



CONCLUSIONS 

We have outlined a numerical approach, canonical di- 
agonalization, for treating a variety of quantum many 
body problems. CD is quite different from most ex- 
isting methods for treating such problems: it does not 
utilize approximate wavefunctions, scmiclassical approx- 
imations, path integrals, perturbation theory, or Monte 
Carlo. Instead, the second quantized Hamiltonian is 
transformed directly, using canonical transformations, to 
put it into a diagonal form. 

We have demonstrated CD on ab initio quantum chem- 
ical calculations for a small molecule. CD appears to 
be quite competitive with the best alternative quantum 
chemical methods, such as the coupled cluster method, 
even in this early stage of its development. Unlike many 
other approaches, CD can be used to treat systems where 
the ground state has a small overlap with the Hartree 
Fock state. It can also be used to remove high energy or- 
bitals from the problem, leaving a smaller problem which 
can be treated with other methods, such as DMRG. Al- 
though we have not yet tested the ability of CD to obtain 
excited states, there is no fundamental difference between 
the ground state and an excited state within CD, and we 
have outlined specific methods to obtain excited states. 

One of the principle future uses of CD could be to 
derive simple model Hamiltonians, much studied in con- 
densed matter physics, directly from ab initio calcula- 
tions. Currently, deriving model Hamiltonians is an art 
which involves educated guesses for the proper model 
forms coupled with the matching of completely separate 
solutions for the ab initio and model systems. CD may 
be able to unify this approach into a controlled single 
procedure. 
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